function outp = ReentryDynamic(state, ctrl)
%% dimensionless [t,rad,lon,lat,V,fpa,azi];
global m S R0 g0 Vs rot
rad  = state(1);
lat = state(3);
V = state(4);
fpa = state(5);
azi = state(6);
alt = (rad-1)*R0;
[rho, sonic] = getAtmosphere(alt);

alpha = ctrl(1);
absBank = ctrl(2);
signOfBank = ctrl(3);
bank = absBank*signOfBank;
Ma = (V*Vs)/sonic;
[CL, CD] = getCLCD(alpha, Ma);
q = 0.5*rho*(V*Vs)^2;
D = q*S*CD/m/g0;
L = q*S*CL/m/g0;

raddot   = V*sin(fpa);
londot   = V*cos(fpa)*sin(azi)/(rad*cos(lat));
latdot   = V*cos(fpa)*cos(azi)/rad;
speeddot = -D - sin(fpa)/(rad^2) + rot^2*rad*cos(lat)*(sin(fpa)*cos(lat)-cos(fpa)*sin(lat)*cos(azi));
fpadot   = L*cos(bank)/V + cos(fpa)*(V/rad - 1/V/rad/rad) + 2*rot*cos(lat)*sin(azi)...
      + rot^2*rad*cos(lat)*(cos(fpa)*cos(lat) + sin(fpa)*cos(azi)*sin(lat))/V;
azidot   = L*sin(bank)/(V*cos(fpa)) + V*cos(fpa)*sin(azi)*tan(lat)/rad +...
    2*rot*(sin(lat)-tan(fpa)*cos(lat)*cos(azi)) + rot^2*rad*sin(azi)*sin(lat)*cos(lat)/(V*cos(fpa));
outp  = [raddot, londot, latdot, speeddot, fpadot, azidot];
end